Critical temperature of a trapped Bose gas: comparison of theory and experiment 
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We apply the Projected Gross-Pitaevskii equation (PGPE) formalism to the experimental problem of the shift 
in critical temperature T c of a harmonically confined Bose gas as reported in Gerbier et al. [Phys. Rev. Lett. 
92, 030405 (2004)]. The PGPE method includes critical fluctuations and we find the results differ from various 
mean-field theories, and are in best agreement with experimental data. To unequivocally observe beyond mean- 
field effects, however, the experimental precision must either improve by an order of magnitude, or consider 
more strongly interacting systems. This is the first application of a classical field method to make quantitative 
comparison with experiment. 
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The shift in critical temperature T c with interaction strength 
for the homogeneous Bose gas has been the subject of numer- 
ous studies spanning almost fifty years since the first calcula- 
tions of Lee and Yang QLQ]. While there is a finite shift to the 
chemical potential in mean-field (MF) theory, the shift of the 
critical temperature is zero id. The leading order effect is due 
to long-wavelength critical fluctuations and is inherently non- 
perturbative. Using effective field theory it was determined 
that the shift is AT C /T® = can 1 ' 3 where n is the particle 
number density, a is the s-wave scattering length, and c is a 
constant of order unity 0]. Until recently results for the value 
of c disagreed by an order of magnitude and even sign, as sum- 
marised in Fig. 1 of JH]. However, two calculations performed 
using lattice Monte Carlo have settled the matter, and confirm 
that the shift is in the positive direction with combined esti- 
mate of c ~ 1.31 ± 0.02 pi [a]. A number of recent improved 
results broadly agree, and useful discussions are provided by 
Andersen 0] and Holzmann et al. jH]. 

The situation for the harmonically confined Bose gas is 
somewhat different. The ideal gas transition temperature and 
de Broglie wavelength at T c are 



mil 



( N 



1/3 



k B ' 



A = 



2irh z 



1/2 



(1) 



VC(3) 

with Co = (tUxUiyUz) 1 / 3 . There is a shift in T c due to finite 
size effects U given by AT C /T C ° ~ -Q.HN~ 1 / 3 Cd/uj with 
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experimentally relevant parameters. The first-order shift in T c 
that survives in the thermodynamic limit is due to mean-field 
effects and has been estimated analytically |flol] . Repulsive 
interactions reduce T c , which can be intuitively understood 
due to a lowering of the peak density of the gas. Next-order 
effects due to fluctuations have been estimated in 
and in general predict an increase in T c from the first order 
result. For a sufficiently wide trap Ref. 11211 estimates 
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with ci = -3.426, c' 2 = -45.86, c 2 ' = -155.0, which for 
a/A < 0.032 predicts a positive shift due to fluctuations. 
The first term is the MF result of ifloll . Recently Zobay and 
co-authors have investigated power law traps with the goal of 
understanding how T c behaves in a smooth transition from 
harmonic trapping to the homogeneous situation 

For a typical BEC experiment, the critical temperature de- 
viates from the ideal gas result only by a few percent. As 
thermometry of Bose gases at this level of accuracy can be dif- 
ficult Oil, until recently the only experimental measurement 
was reported by Ensher et al. with AT C /T C ° = -0.06 ± 0.05 
lfl7ll . However, in 2004 the Orsay group reported precise mea- 
surements of the critical temperature for a range of atom num- 
bers, and compared their results to the first-order MF estimate 
of Poll . While in agreement, the theoretical results lie near the 
upper range of the experimental error bars. 

Previously one of us used the classical field projected 
Gross-Pitaevskii equation (PGPE) formalism OH 111 iSl to 
give an estimate of the shift in T c of the homogeneous Bose 
gas I2TI1 . which was found to be in agreement with the Monte 
Carlo calculations H, @1, The PGPE is a dynamical non- 
perturbative method, with the only approximation being that 
the highly occupied modes ((Nk) 1) of the quantum Bose 
field are well-approximated by a classical field evolved ac- 
cording to the GPE. Related classical field approaches have 
been considered by a number of authors, including Kagan and 
co-workers i22ll . Sinatra et al. i23l 12411 . Rzazewski and co- 
workers ll25l l26ll . 

Here we use an extension of the PGPE for harmonically 
trapped gases l27ll to calculate the shift in T c for the experi- 
ment of Gerbier et al, and in particular focus on the compet- 
ing effects of mean-field and critical fluctuations. The PGPE 
in dimensionless units is 
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(2) where ^> is the classical field, V = (A 2 x 2 + A 2 ?/ 
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The nonlinearity is C n \ = Ni,Uo/Tiu) z XQ 
^jTi/2mu> z and r = u z t. For 
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and X x ,y — (*Jx,y/^2 
where the unit of length is xq 
the harmonic trap the Bose field is expanded on a basis of 
harmonic oscillator eigenstates, with the cutoff energy E cut 
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FIG. 1: (color online) Comparison of condensate fraction versus 
temperature normalised to T c for the PGPE method and the HFB- 
Popov calculation for = V°, Xy = 1 and E cut = Slhui,. PGPE 
method: C n i = 500 (crosses), C n \ = 2000 (circles), C n i = 10000 
(stars). Dashed line: HFB-Popov results. Solid line: exact result for 
Cni = 0. Inset: AT = Tpgpe — Thfb at fixed condensate fraction 
for C n i = 2000 and 10000, indicating a maximum shift near T c . 



determined by the occupation number condition. The projec- 
tion operator V{F} projects the function F onto the harmonic 
oscillator modes with energy less than E cut . 

The dynamical PGPE system represents a microcanonical 
ensemble, and will evolve any random initial conditions to 
thermal equilibrium defined by the integrals of motion lfl9tl . 
For a cylindrically symmetric harmonic trap these are the to- 
tal number of particles, the energy, and the component of the 
angular momentum along the symmetry axis. Once in equi- 
librium, we use the assumption of ergodicity to accurately de- 
termine the condensate fraction 12711 . and the temperature T 
and chemical potential /i& |2lll . By varying the initial state 
energy we measure the dependence of condensate fraction on 
temperature. 

As an initial investigation into critical fluctuations, in Fig.Q] 
the results of the PGPE calculations from 1I27I1 are compared 
with a self-consistent mean-field calculation in the Popov ap- 
proximation to the Hartree-Fock-Bogoliubov (HFB) formal- 
ism (see e.g. l28ll ). In order to make a direct comparison, the 
HFB-Popov calculation is performed in the same basis as the 
dynamical PGPE calculations, and we use the equipartition 
distribution Nk = fc^T / (e*. — fJ.) for the quasi-particle occu- 
pations. (This is the high temperature limit fcgT 3> ek of the 
Bose-Einstein distribution applicable to classical fields). For 
smaller values of C D \ the HFB-Popov theory agrees with the 
classical field calculation, however for larger values there is a 
distinct difference which we attribute to critical fluctuations. 
We have repeated these calculations using gapless implemen- 
tations of HFB theory j29[| and found that they are little differ- 
ent from the results calculated using HFB-Popov. Our results 
demonstrate that critical fluctuations have a measurable effect 
for the PGPE system. However this is an idealised calculation 
— to be quantitative we must make a connection between the 
PGPE method and the recent experiment l30ll . 
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TABLE I: Input parameters for the PGPE simulations. The chemical 
potential m and the shift of the cutoff energy S c are outputs parame- 
ters measured at the critical point. 



Gerbier et al. IB0H trap 8 Rb atoms in a cylindrically sym- 
metric harmonic potential with (u> x>y ,u) z ) = 2ir x (413, 8.69) 
Hz giving X XtV = 47.52. For total numbers of atoms iV t ot 
ranging from 2.5 x 10 5 to 2.5 x 10 6 , the critical point was 
identified by reducing the final rf frequency of the evaporative 
cooling, and identifying the point that the condensate fraction 
became measurable (see Fig. 2 of l30ll .) We perform numeri- 
cal simulations in a similar manner. We choose relevant simu- 
lation parameters and dynamically evolve the system to equi- 
librium for a range of energies. We identify the critical point 
from the number of condensate particles and determine the 
number of particles above the cutoff using a self-consistent 
semi-classical approximation for the high-energy modes as 
described below. This gives us a set of points (N c , T c ) to be 
compared with the experimental data. 

To simulate the experiments of Gerbier et al. using the 
PGPE we need to choose both an energy cutoff E cut and a 
number of particles below the cutoff Nb to simulate so that 
the occupation number condition is satisfied. However, any 
final result should be insensitive to the exact value of the 
cutoff that is chosen. A priori estimates for our simulation 
parameters were determined from the Bose-Einstein distribu- 
tion of quantum orbitals of an ideal trapped gas at the critical 
temperature, and are summarised in Table U For the smaller 
clouds we chose an energy cutoff such that (Nk) > 5. For the 
large clouds this leads to correspondingly larger basis sets that 
become computationally prohibitive, and for these we chose 
{Nk) > 7.5. In principle we could use this occupation con- 
dition for all simulations, however the two calculations at the 
crossover point (iV t ot = 2.5 x 10 6 ) enables us to verify that 
our calculations are insensitive to the exact value of the en- 
ergy cutoff. We use the PGPE to evolve randomised initial 
states to equilibrium and measure the condensate number No, 
chemical potential /it, temperature T, and density n&(x) for 
each set of parameters. 

In Fig. |2 a) we plot the condensate number versus tempera- 
ture for the N® ot = 4x 10 6 data set and find there is no sharp 
transition. This is because we are only considering the atoms 
below the cutoff. As the majority of atoms in the full system 
are above the cutoff and Nq is of order a few hundred particles 
for all the data points on this graph, these simulation results 
all lie close to the critical point. To determine a single critical 
point from each data set we plot on the same graph the corre- 
sponding condensate number for the finite-sized ideal gas at 
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FIG. 2: (color online) Determination of the critical number and tem- 
perature for the simulation with AT t ° ot = 4.0 x 10 6 . (a) iVo vs T 
for: classical field (crosses, plusses, solid line), HFB-Popov (circles, 
dashed line). The number of condensate atoms for the ideal gas at the 
critical temperature is the dot-dash line, (b) Total number of atoms 
for: classical field (crosses, plusses, solid line), HFB-Popov (circles, 
dashed line). Critical number versus temperature for the finite-sized 
ideal gas is the dot-dash line. For both (a) and (b) the solid and 
dashed lines are polynomial fits to the data. 



the same critical temperature. We choose the intersection of 
these two curves to identify the critical point, and have veri- 
fied that the occupation number condition is satisfied here. 

To relate these results back to the full experimental system 
we assume that the classical field and the above cutoff thermal 
cloud are weakly-coupled systems in equilibrium, with the 
same temperature and chemical potential. The thermal cloud 
exists in the potential of the trap plus time-averaged classical 
field density n& (x) that is determined from the PGPE simula- 
tions. To solve for the above cutoff thermal cloud we make use 
of the self-consistent Hartree-Fock approximation, which pro- 
vides an accurate description of the modes above E cut . The 
above cutoff density is determined by the self-consistent solu- 
tion of 
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£W(p,x) = p 2 /2m + T4ra P (x) + 2[/ K(x) + n„(x)], (5) 

where £"hf(p,x) is the Hartree-Fock energy. In this pro- 
cedure the contribution of the above cutoff density n a (x) 
to the effective potential for the classical field is neglected. 
This is justified as we find that near the critical point n a (x) 
is approximately flat in the region where the 7i&(x) is non- 
zero. However, the uniform energy shift of this interaction 
must be included in the chemical potential used in Eq. © as 
(J, = fib + 2n (0)E/o- Another important correction accounts 
for the shift in the energy of the highest oscillator modes in 
the classical field from E cut due to interaction effects so that 
the integral in Eq. is over the correct region of phase space. 
We do this by assuming that the highest energy modes of the 
classical field are single-particle in nature, and are shifted by 
a constant amount 5 C . We fit the time-averaged occupation of 
these modes to (iVjt) = T j (ejj + 5 C — fib)- The lower limit of 
the integral in Eq. is then Eq = E cxlt + S c + 2n a (0)Uo to 
account for the mean-field of the thermal cloud. 
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FIG. 3: (color online) Comparison of theoretical calculations with 
experiment. The main figure plots T c vs N c , whereas the inset plots 
the shift of T c against the relevant small parameter a/Ao. Experi- 
mental results: data (open circles), one a fit (grey area). Theoret- 
ical results for T c : ideal gas (dot-dash), Al (dotted), A2 (dashed), 
MF-GPE (crosses), MF-HFBP (dots), PGPE (plusses). Solid lines 
through the data points are polynomial fits. Al is not shown in the 
main figure for clarity. 



We have also calculated T c using other methods for com- 
parison, as summarised below: 

1. Al: This is the first order analytic estimate of Giorgini et 
al. iflOn . which is the first term of Eq. ©. 

2. A2: This is the full second order result of Eq. (|2]). How- 
ever, the validity condition for this result (Eq. (7.2) of lfl~2tl ) 
requires the trap to be "sufficiently broad", and this is strongly 
violated for this experiment. This essentially says that the 
semi-classical approximation is not valid for the lowest en- 
ergy modes of this strongly elongated system. 

3. MF-GPE: The GPE is solved numerically using a varia- 
tional Gaussian ansatz, and the thermal cloud calculated us- 
ing a semi-classical approximation ifioll . At each tempera- 
ture the condensate and non-condensate are determined self- 
consistently with a fixed number of particles, and the criti- 
cal temperature is where the condensate fraction decreases to 
zero. 

4. MF-HFBP: We fix the condensate fraction, and determine 
the temperature that gives an appropriate self-consistent con- 
densate mode and thermal density. We have verified the re- 
sults are unchanged for equipartition or Bose-Einstein statis- 
tics. We use the same procedure as for the PGPE calculation 
to determine the critical point, the above cutoff density and 
the total atom number. An illustrative set of data is displayed 
in Fig. |2] 

In Fig. [3] we compare these theoretical results with the 
PGPE and experimental data. The MF Al estimate was shown 
in i30ll and is within the experimental error bars. However, our 
more accurate MF-GPE calculation gives a greater value of T c 
at larger atom numbers, agreeing with the mean-field results 
of Houbiers et al. ifTHl . However, the MF-HFBP result, which 
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presumably is an even better mean-field calculation, is quite 
different and towards the lower end of the experimental error 
estimate. 

The predicted effect of critical fluctuations il~i1,[l2tl is to fur- 
ther increase T c . The non-perturbative A2 estimate lies at the 
boundary of experimental error, but as mentioned earlier this 
result does not satisfy the validity requirement for this exper- 
iment. The PGPE calculation, which includes all the physics 
of the MF-HFBP calculation as well as critical fluctuations, is 
measurably different. Arguably it is in best agreement with 
the experimental data. However, both the PGPE and MF- 
HFBP calculations lie within the error bars, suggesting that 
experimental precision must improve by an order of magni- 
tude in order to distinguish these predictions. 

The inset of Fig. [3] shows the PGPE shift as a function of 
a/Ao and in comparison with the results of Eq. (O and the 
experimental data. The second order term is almost constant 
over the experimental range of a/Ao and so cannot distinguish 
the presence or otherwise of any logarithmic term. We note 
that the finite-size shift is subtracted from the PGPE and ex- 
perimental data for this comparison. 

We have also translated data for parameters as in Fig. Q] 
but with C n i = 5000 to realistic experimental values, and 
found that for 10 7 atoms of 87 Rb in a TOP trap with a 40 Hz 



radial frequency that the difference between the MF-HFBP 
and PGPE results is of order 3%. Thus we suggest that for 
currently accessible experimental conditions it will be neces- 
sary to either make use of Feshbach resonances to probe more 
strongly interacting regimes, or to move to traps flatter than 
harmonic to be able to distinguish these theories in the lab. 

In conclusion we have performed a careful theoretical anal- 
ysis of the experiment on the shift in critical temperature of 
a trapped Bose gas reported in Gerbier et al. IT30ll . We have 
determined that earlier calculations based on mean-field the- 
ory and the local-density approximation are inappropriate for 
this experiment, and make predications for T c outside the ex- 
perimental error bars at larger atom numbers. We have ap- 
plied non-perturbative classical field theory to this problem, 
and described how to incorporate the physics of the above- 
cutoff atoms in equilibrium. The results include the effect of 
critical fluctuations, and give the best agreement with exper- 
imental observations. Our results indicate the precision re- 
quirements for experiments to investigate beyond mean-field 
effects on the critical temperature. 
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